#include "bi_calcBlocksize.hpp"
#include "cacheInfo.hpp"

#include <cmath>
#include <limits>

namespace hms::swe
{

namespace block {

double getApproxSize(){
/* the cache commands currently rely on grep and lscpu, or sysconf */
#ifdef HMS_IS_POSIX
	/* The assumption is that the optimal number of cells per block is
	 * proportional to the last level cache size available to each thread. */
	Index cpt { hms::cachePerThread() };
	/* All cell-based variables need to be stored for each cell. This factor
	 * is empirical, and subject to change, if the solver uses more or fewer
	 * variables per cell */
	static constexpr int dataSizePerCell { 16 };
	/* Assuming a block shape close to a square, the number of cells per block
	 * is approximately the square of the block size. */
	return { std::sqrt(cpt) / dataSizePerCell };
#else 
	return 64;
#endif
}

Index calculateSize_threshold(
	Index meshsize_dim,
	Index approxSize,
	double threshold
){
	Index
		minSize { static_cast<Index>( std::ceil(threshold * approxSize) ) },
		tweakedSize { static_cast<Index>(approxSize) },
		maxRem { std::numeric_limits<Index>::min() },
		rem;

	for ( Index i {static_cast<Index>(approxSize)}; i>=minSize; --i ){
		rem = meshsize_dim % i;
		if ( rem == 0 ){
			tweakedSize = i;
			break;
		}
		if ( rem > maxRem ){
			maxRem = rem;
			tweakedSize = i;
		}
	}
	return tweakedSize;
}

Index calculateSize_nextNeighbour(
	Index meshsize_dim, double approxSize
){
	Index
		nBlocks { static_cast<Index>(
			std::ceil( meshsize_dim / approxSize )
		)},
		size_ceil { static_cast<Index>( std::ceil(  meshsize_dim / nBlocks ) )},
		size_floor{ static_cast<Index>( std::floor( meshsize_dim / nBlocks ) )},
		rem_ceil  { meshsize_dim % size_ceil  },
		rem_floor { meshsize_dim % size_floor };

	if ( rem_ceil == 0 )
		return size_ceil;
	else if ( rem_floor == 0 )
		return size_floor;
	else if ( rem_ceil >= rem_floor )
		return size_ceil;
	else
		return size_floor;
}

} // namespace block

Index calculateBlocksize(
	Index meshsize_dim, RemAdjust alg, double threshold
){
	double approxSize { block::getApproxSize() };
	/* Take the remainder size into account.
	 * There are two competing optimisation targets:
	 * 1. The block size should be as close to approxSize as possible
	 * 2. The remainder size should be zero or as close to approxSize
	 *    as possible.
	 * As a condition, the block size must never exceed approxSize. */
	Index tweakedSize { [&](){
		using R = RemAdjust;
		switch (alg){
			case R::threshold:
				return block::calculateSize_threshold(
					meshsize_dim, approxSize, threshold
				);
			case R::nextNeighbour:
				return block::calculateSize_nextNeighbour(
					meshsize_dim, approxSize
				);
			case R::none: [[fallthrough]];
			default:
				return static_cast<Index>(approxSize);
		}
	}() };
	/* Account for eventualities and define a minimum and maximum size.
	 * Currently (2021-11-25), the largest cache size on any cpu is 144 MB. */
	static constexpr Index
		minSize { 32 }, maxSize { 768 };
	return std::min( std::max(minSize, tweakedSize), maxSize );
}

} // namespace hms::swe
